include("TKTD_04_boilerplate_dl_ni.jl"); # set up simulation and load default parameters for d magna

  Activating project at `c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations`

theme(:ggplot2) # ggplot theme works better in jupyter 
default(leg=false)

# using the longispina AMP parameters as default
default_flea_params = copy(default_fleas[:DL_AMP])
# turn off ageing and adaptive plasticity by default
default_flea_params[:h_a_accel] = 0.
default_flea_params[:e_50_adpt] = 1e-10
default_flea_params[:beta_e_adpt] = 10.

# phytoplankton (food) parameters
phyto_params = [copy(default_phyto.Rsub)]

1-element Vector{OrderedDict{Symbol, Any}}:
 OrderedDict(:species => "Rsub", :mu_max => 1.78, :q_max => 0.018, :q_min => 0.0011, :q_init => 1.0, :v_max => 0.062, :k_s => 0.0625, :m_max => 0.0, :I_opt => 120.0, :T_opt => 28.0…)

Implementation verification

The following code has the purpose to validate the DEB-IBM implementation for each submodel.
An abbreviated version will be provided in the TRACE documentation.

To facilitate comparison of model output with analytical solutions, some changes to the simulation settings are made:

default_global_params[:A_inoculate] = [100.] # => f=1
default_global_params[:renew_medium] = Int64[] # no renewals
default_global_params[:remove_juveniles] = 0:7 |> collect; # remove juveniles immediately
default_global_params[:age_init_max] = 1e-7; # start directly at birth

default_global_params[:data_recording_interval] = 1.;

DEB

Embryonal development time ("age at birth")

According To Kooijman (2010), age at birth (i.e., embryonal development time) can be predicted analytically when $g$ is small and $\dot{k}_M$ is large:

$$ a_b = \frac{3l_b}{\dot{k}_M g}\left( 1 + \frac{g}{f} \right) $$

, with $L_{max}=\kappa \frac{\left\{ \dot{p}_{A,max} \right\}}{\left[ \dot{p}_M \right]}$, $g=\frac{ [E_G]\dot{v} }{ \kappa \left\{ \dot{p}_{A,max} \right\} }$ and $\dot{k}_M=\frac{\left[ \dot{p}_M \right]}{E_G}$.

This equation also requires that $f$ is equal to scaled reserve density at birth, $e_b$.
$l$ is length $L$ scaled by maximum length $L_{max}$.
We know all the inputs except for $l_b$. So we cannot calculate $a_b$ directly, but we can check for consistency by inserting $l_b$ from the simulation output table, calculating $a_b$ correspondingly, and comparing with $a_b$ in the simulation output table (which is the actual age at birth that results from simulating embryonal development).

function L_max(
        kappa::Float64,
        p_A_m::Float64,
        p_M_V::Float64
    )
    return kappa*(p_A_m/p_M_V)
end

function g(
        E_G::Float64, 
        v::Float64, 
        kappa::Float64, 
        p_A_m::Float64)
    return (E_G*v)/(kappa*p_A_m)
end

function k_M(
        p_M_V::Float64,
        E_G::Float64
    )
    return p_M_V/E_G
end

k_M (generic function with 1 method)

"""
Calculate age at birth for small g and large k_M. 
$(TYPEDSIGNATURES)
"""
function age_at_birth(
        L_b::Float64,
        L_max::Float64,
        k_M::Float64,
        g::Float64,
        f::Float64
    )
    let l_b=L_b/L_max
        return ((3l_b)/(k_M*g))*(1+(g/f))
    end
end

age_at_birth

The code below will run simulation of a life-table experiment for a single test animal, and show output for length and age at birth for a single individual.

flea_params = deepcopy(default_flea_params) 

out = life_table(default_global_params, phyto_params, flea_params; n_reps=1)
L_b_sim, a_b_sim = @subset(out.repro, :t_day.==minimum(:t_day))[1,[:L, :age]]

DataFrameRow (2 columns)

	L	age
	Any	Any
1	0.017281	3.25

This code cell will show the analytical solution, making use of the simulated $L_b$. The value that is printed should be equal to "age" in the table above.

L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], Float64(flea_params[:p_M_V_0]))
g_anl = g(Float64(flea_params[:E_G_0]), flea_params[:v], flea_params[:kappa_0], flea_params[:p_A_m_0])
k_M_anl = k_M(Float64(flea_params[:p_M_V_0]), Float64(flea_params[:E_G_0]))
a_b_anl = age_at_birth(L_b_sim, L_max_anl, k_M_anl, g_anl, 1.)

3.2678204441392853

The values are fairly close, but it might be useful to generate repeated predictions and also use different input values for $\kappa$ and $\dot{v}$ to get a better idea.

"""
Get simulated and analytically calculated age at bith, given values for kappa and v.
$(TYPEDSIGNATURES)
"""
function evaluate_a_b(kappa::Float64, v::Float64)
    flea_params = deepcopy(default_flea_params) 
    # smallest E_G_0 value that doesn't break the model (k_M becomes large as a result of this)
    flea_params[:kappa_0] = kappa
    flea_params[:v] = v

    out = life_table(default_global_params, phyto_params, flea_params; n_reps=1)
    L_b_sim, a_b_sim = @subset(out.repro, :t_day.==minimum(:t_day))[1,[:L, :age]]

    L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], flea_params[:p_M_V_0])
    g_anl = g(flea_params[:E_G_0], flea_params[:v], flea_params[:kappa_0], flea_params[:p_A_m_0])
    k_M_anl = k_M(flea_params[:p_M_V_0], flea_params[:E_G_0])
    a_b_anl = age_at_birth(L_b_sim, L_max_anl, k_M_anl, g_anl, 1.)
    return a_b_sim, a_b_anl
end

vec_a_b_sim = []
vec_a_b_anl = []
kappa_sim = []
for kappa in range(0.1, 0.9, length=10)
    try
        a_b_sim, a_b_anl = evaluate_a_b(kappa, default_flea_params[:v])
        push!(vec_a_b_sim, a_b_sim)
        push!(vec_a_b_anl, a_b_anl)
        push!(kappa_sim, kappa)
    catch
        nothing
    end
end
p = scatter(
    kappa_sim, vec_a_b_anl, 
    xlabel=L"\kappa\ (-)", ylabel="Embryonal \n development time (d)", 
    leg=:outertopright, label="Analytical solution", 
    size=(800,250), leftmargin=5mm, bottommargin=5mm
)
scatter!(kappa_sim, vec_a_b_sim, label="Simulation output")
savefig(plot(p, dpi=300), "figures/implementation_verification_embryonal_development_time.png")
p

vec_a_b_sim = []
vec_a_b_anl = []
v_sim = []

for v in range(0.0022, 0.02, length=10)
    try
        a_b_sim, a_b_anl = evaluate_a_b(default_flea_params[:kappa_0], v)
        push!(vec_a_b_sim, a_b_sim)
        push!(vec_a_b_anl, a_b_anl)
        push!(v_sim, v)
    catch
        nothing
    end
end
scatter(
    v_sim, vec_a_b_anl, 
    xlabel=L"\dot{v}\ (cm\ d^{-1})", ylabel="Age at birth (d)", 
    leg=:outertopright, label="Analytical solution", 
    size=(800,250), leftmargin=5mm, bottommargin=5mm
)
scatter!(v_sim, vec_a_b_sim, label="Simulation output")

The simulation output follows the same pattern as the analytical solution. There is a consistent and very marginal difference in the absolute values of $a_b$. It is not entrirely clear what causes it, but the differences are probably small enough to be irrelevant, and occur for parameter combinations which are biologically irrelevant.\ Note that the analytical solution is only valid under specific assumptions. When other (potentially realistic) input parameters are used, it is not guaranteed that the same level of agreement will be reached.

Maximum length

Maximum length (at constant and abundant food) can be directly calculated from parameters: $$ L_{max} = \kappa \frac{ \left\{ \dot{p}_{A,max} \right\} }{ \left[ \dot{p}_M \right] } $$

To verify that $L_{max}$ is replicated correctly, we need to set assimilation efficiency constant and calculate it internally from maximum ingestion and assimilation rates.

flea_params = deepcopy(default_flea_params) 

flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
# increase t_max to be sure that the analytical L_max will be reached
global_params[:t_max] = 56;
# calculate L max analytically
L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], flea_params[:p_M_V_0])

0.054705986990994986

out = life_table(global_params, phyto_params, flea_params; n_reps=1)

p = @df out.growth plot(
    :t_day, :L, 
    xlabel="Time since birth (d)", ylabel="Structural length (cm)", 
    leg=:bottomright, label="Simulation output",
    ylim=(0, 0.06)
    )
hline!([L_max_anl], label="Analytically predicted maximum length")
savefig(plot(p, dpi=300), "figures/implementation_verification_maximum_length.png")
p

Feeding

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_50_adpt] = 0.75
flea_params[:beta_e_adpt] = 2.0
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.1]
#global_params[:renew_medium] = [0, 2, 5]

1-element Vector{Float64}:
 0.1

out = life_table(global_params, phyto_params, flea_params; n_reps=1);

We simulated food being added daily, but never being removed by medium renewals. Since a single animal cannot ingest that much, food accumulates in the medium:

@df out.repro plot(
    plot(:t_day, :X_p , linetype=:steppost), size=(400, 200), xlabel="Time (d)", ylabel="Food density (J/L)",
    xticks=0:7:21
    )

We can use this output to check the shape of the scaled functional response, which should perfectly follow the Type II functional response equation.
This can be compared with the functional response plotted on top:

# calulating the functional response directly
xvals = range(extrema(out.repro.X_p)..., length=100)
f_anl = @. xvals / ((flea_params[:J_X_A_m_0]/flea_params[:F_m_0])+xvals)
# comparison with simulation output
   
plot(xvals, f_anl, xlabel="Food density (J/L)", ylabel="f", leg=:bottomright, label="Analytical solution", size=(400,200))
@df out.repro scatter!(:X_p, :f, label="Simulation output")

This shows that the functional response is implemented correctly, but it might be useful to check that the scaling by surface area is also correct. We can effectively fix $f$ to a value of 1 by using a very high food density, then plot ingestion rate against surface area and should see a perfectly linear relationship:

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [1000.]
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro scatter(:L.^2, :J_X, xlabel=L"L^2 (cm^2)", ylabel=L"\dot{J_X} (J/d)", size=(400,300))

Lastly, we need to check that the external food density is updated correctly. We do so by calculating the differential along the vector of output food densities and comparing this with the output ingestion rate, which, in the absence of medium renewals, should be mirror images.
They might not be perfectly reciprocal because neonates are only removed daily and will also ingest some food, but the difference should not exceed a few percent of the total change in food density.
For this, we turn off daily addition of food, so food densities should decline continuously.

flea_params = deepcopy(default_flea_params) 
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1/24
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
out = life_table(global_params, phyto_params, flea_params; n_reps=1);

┌ Warning: Inference of survival will currently fail for data_recording_interval !=1.
└ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:97

dX = diffvec(out.repro.X_p) # change in external food concentration
JX = out.repro.J_X # ingestion rate
ratio = abs.(dX ./ JX) # ratio between both
extrema(ratio[2:end]) # show extrema of ratios, excluding first timepoint (=always 0)

(1.0239657417767893, 1.04331586949277)

Adaptive plasticity

For adaptive plasticity, there are no analytical solutions that would be useful for implementation verification.
We do a quick check that $\kappa$ increases with increasing $e$ and approaches $kappa_0$.

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_50_adpt] = 0.5
flea_params[:beta_e_adpt] = 4.0
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.1]
#global_params[:renew_medium] = [0, 2, 5]

1-element Vector{Float64}:
 0.1

out = life_table(global_params, phyto_params, flea_params; n_reps=1);

Note that it is not guaranteed that $\kappa_0$ will be reached. This depends on the value of $\beta_{e,adpt}$.
$kappa_0$ (horizontal line) should however be approached as $e$ increases.

@df out.repro plot(:e, :kappa, xlim=(0,1), ylim=(0,1), label=L"\kappa_e", xlabel="e", ylabel=L"\kappa", marker=true)
hline!([flea_params[:kappa_0]], label=L"\kappa_0")

Reproduction

Reproduction should occur in discrete intervals, dictated by the time between molts. Note that the plot below indicates reproduction at the time-points when neonates hatch (after embryonal development has been completed), not at the time-points of oviposition.

flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro plot(:t_day, :cum_repro, linetype=:steppost, xlabel="Time (d)", ylabel="Cumulative reproduction", leg=false)

┌ Warning: Inference of survival will currently fail for data_recording_interval !=1.
└ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:97

If we increase $\tau_{molt}$, reproduction should occur less frequently. ($t_{max}$ is increased to demonstrate this).
Note that the values used below are not biologically realistic for Daphnia, but only used to demonstrate that increase of the molting period has the expected effect

flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:tau_molt] = 7
flea_params[:h_a] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
global_params[:t_max] = 56
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro plot(:t_day, :cum_repro, linetype=:steppost, xlabel="Time (d)", ylabel="Cumulative reproduction")

┌ Warning: Inference of survival will currently fail for data_recording_interval !=1.
└ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:97

The energy in the reproduction buffer should increase continously during adulthood, and decrease in a discrete manner.

flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
out = life_table(global_params, phyto_params, flea_params; n_reps=1)
@df out.repro plot(:t_day, :E_R, linetype=:steppost, xlabel="Time (d)", ylabel=L"E_R\ (J)")

┌ Warning: Inference of survival will currently fail for data_recording_interval !=1.
└ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:97

Starvation mortality

Survival is simulated at two food levels for three different values of the starvation parameter $\beta_e$.

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56

outs = LifeTableDataset[]
e_starve_sweep = [0.125, 0.25, 0.5]
for e_starve in e_starve_sweep
    for A_0 in [0.1, 0.4]
        global_params[:A_0] = [A_0]
        flea_params[:e_starve] = e_starve
        flea_params[:beta_e] = 0.5
        out = life_table(global_params, phyto_params, flea_params; n_reps=10)
        out.survival[!,:e_starve] .= e_starve
        out.survival.food_level .= A_0
        push!(outs, out)
    end
end
out = concat(outs)
typeof(out)

LifeTableDataset

plot(
    (@df @subset(out.survival, :food_level.==0.4) plot(
        :t_day, :survival, group=:e_starve, linetype=:steppost,
        xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)", 
        legend=false,
        leftmargin=5mm, bottommargin=5mm, xticks=0:14:56, title="High food \n (0.4 mg dm/day)"
        )),
    (@df @subset(out.survival, :food_level.==0.1) plot(
        :t_day, :survival, group=:e_starve, linetype=:steppost,
        xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)", 
        legend_title="Starvation \n parameter (-)", 
        leftmargin=5mm, bottommargin=5mm, xticks=0:14:56, title="Low food \n (0.1 mg dm/day)",
        titlefontsize=12, legend=:outertopright
        )), 
        lw=1.5, background=:transparent, labelfontsize=14, layout=@layout([a{0.4w} b{0.6w}]),
        size=(1200,450), topmargin=5mm, ylim=(0,1)
)

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_starve] = 0.25

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.1]
global_params[:t_max] = 56

outs = LifeTableDataset[]
beta_e_sweep = [0.01, 0.1, 2.0]
for beta_e in beta_e_sweep
    flea_params[:beta_e] = beta_e
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.survival.beta_e .= beta_e
    push!(outs, out)
end
out = concat(outs)
typeof(out)

LifeTableDataset

@df out.survival plot(
    :t_day, :survival, group=:beta_e, 
    legend=:outertopright, linetype=:steppost,
    xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)", 
    legend_title=L"Starvation\ parameter\ (-)", size=(800,300),
    leftmargin=5mm, bottommargin=5mm, xticks=0:14:56
)

Ageing

The ageing submodel is checked for plausability by altering the ageing acceleration $\ddot{h}_a$ and comparing survival rates.

flea_params = deepcopy(default_flea_params) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56

outs = LifeTableDataset[]
h_a_sweep = [0, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1]
for h_a in h_a_sweep
    flea_params[:h_a_accel] = h_a
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.survival.h_a .= h_a
    push!(outs, out)
end
out = concat(outs)
typeof(out)

LifeTableDataset

As we would expect, high values for $\ddot{h}_a$ cause individuals to die from ageing at a youger age.

@df out.survival plot(
    :t_day, :survival, group=:h_a, 
    legend=:outertopright, linetype=:steppost,
    xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)", 
    legend_title=L"Ageing\ acceleration\ (1/d^2)", size=(800,300),
    leftmargin=5mm, bottommargin=5mm, xticks=0:14:56
)

TKTD

Toxicokinetics

TK is controlled by a single parameter $\dot{k}_d$. The TK module is checked for plausibility by turning off TD (setting all parameters to 0.) and only simulating TK.
Comparison with analytical solutions is possible, but currently considered too time-intesive for the purpose of validating this implementation.

The plot below shows that, for a $k_d$ value of 1 $d^{-1}$, the substance is taken up very fast.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56.

flea_params = deepcopy(default_fleas.DM_AMP) 
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:TKTD][1].TK[1].params = [1.]

outs = LifeTableDataset[]
C_sweep = [0., 10., 100.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, linetype=:steppost, leg=:outertopright, 
    xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
    size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5
)

Lowering $k_d$ to 0.1 leads, expectedly, to a more gradual uptake.

flea_params[:TKTD][1].TK[1].params = [.1]

outs = LifeTableDataset[]
C_sweep = [0., 10., 100.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, linetype=:steppost, leg=:outertopright, 
    xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
    size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5
)

flea_params[:TKTD][1].TK[1].params = [.05]
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, leg=:outertopright, 
    xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
    size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5, linetype=:steppost
)

outs = LifeTableDataset[]
kd_sweep = 10 .^ range(log10(0.005), log10.(1), length=100)
for k_d in kd_sweep
    global_params[:C][1] = 100.
    flea_params[:TKTD][1].TK[1].params = [k_d]
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.k_d .= k_d
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df @subset(out.repro, round.(:t_day, sigdigits=2).==21) plot(:k_d, [c[1] for c in :C_int_1], xscale=:log10, ylim=(0,101),
    xlabel=L"\dot{k}_d\ (d^{-1})", ylabel="Concentration", label="Scaled damage at day 21"
)
hline!([100.], leg=:bottomright, label="External concentration")

Linking/unlinking of TK models

TK models can be linked for different PMoAs, or independent of each other.
When link_TK=true, the damage values for all PMoAs should be identical.

flea_params[:TKTD][1].TK[1].params = [.05]
flea_params[:TKTD][1].link_TK = true
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
# matrix of damage values per TK component
damage_matrix = vcat(out.repro.C_int_1...)
# this statement needs to be true
damage_matrix[:,1]==damage_matrix[:,2]==damage_matrix[:,3]==damage_matrix[:,4]==damage_matrix[:,5]

true

When link=false, damage dynamics are different for every PMoA.

flea_params[:TKTD][CU].TK[1].params = [.05]
flea_params[:TKTD][CU].TK[2].params = [1.0]
flea_params[:TKTD][CU].TK[3].params = [0.]
flea_params[:TKTD][CU].link_TK = false
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day);

The figure below should show slow TK for TK component 1, fast TK for TK component 2, and no increase in damage for TK component 3.

@df out.repro plot(
    groupedlineplot(:t_day, [c[1] for c in :C_int_1], :C, marker=true, title="TK component 1", ylim=(0,100)),
    groupedlineplot(:t_day, [c[2] for c in :C_int_1], :C, marker=true, title="TK component 2"),
    groupedlineplot(:t_day, [c[3] for c in :C_int_1], :C, marker=true, title="TK component 3", ylim=(0,100)),
    xlabel="Time (d)", ylabel="Scaled \n damage", bottommargin=5mm, leftmargin=5mm, size=(1000,600),
    background_color=:transparent
)

Toxicodynamics

For the TD module, we first check for the presence/absence of effects on reproduction, growth and survival for each mode of action. Reproduction should be affect by all modes of action. Growth should be affected by all by reproduction efficiency. Maintenance costs and Assimilation efficiency should affect survival by triggering starvation mortality

Growth costs

Top row shows carapace length (cm) over time, bottom row cumulative reproduction. Each column is an exposure concentration.

outs = LifeTableDataset[]
for C in C_sweep
    global_params[:C][CU] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro[!,:C] .= C
    out.growth[!,:C] .= C
    out.survival[!,:C] .= C
    push!(outs, out)
end

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.
flea_params[:TKTD][CU].TD = default_flea_params[:TKTD][1].TD[[G]]
flea_params[:TKTD][CU].TD[1].DRC.params = (100., 2., 5.)
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
    global_params[:C][CU] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro[!,:C] .= C
    out.growth[!,:C] .= C
    out.survival[!,:C] .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

combine(groupby(out.growth, :C)) do df
    L0 = df.carapace_length[1]
    L7 = df[argmin(abs.(df.t_day .- 7)),:carapace_length]
    return DataFrame(
        SGR = (log(L7)-log(L0))/7
    )
end |> x-> plot(x.C, x.SGR, xlabel="External concentration", ylabel="7-day somatic growth rate", xscale=:log10, lw=2, title="Effects on somatic growth rate", marker=true)

Maintenance costs

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.0
flea_params[:TKTD][CU].TD = default_flea_params[:TKTD][1].TD[[M]]
flea_params[:TKTD][CU].TD[1].DRC.params = (100., 2., 5.)
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true

outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

combine(groupby(out.growth, :C)) do df
    L0 = df.carapace_length[1]
    L7 = df[argmin(abs.(df.t_day .- 7)),:carapace_length]
    return DataFrame(
        SGR = (log(L7)-log(L0))/7
    )
end |> x-> plot(x.C, x.SGR, xlabel="External concentration", ylabel="7-day somatic growth rate", xscale=:log10, lw=2, title="Effects on somatic growth rate", marker=true)

Assimilation efficiency

Assimilation efficiency causes effects on all variables, including survival if $\dot{h}_{s,e}$ is large enough.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true

flea_params[:TKTD][CU].TD[A].DRC.params = (100., 2., 1.)

outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

Reproduction efficiency

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[R].DRC.params = (100., 2., 1.)

outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

Same for slow TK to check that the connection between TK and TD works as desired.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.01]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
    LL3GM,
    (100., 2., 1.)
)

outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

Somatic investment ratio

Increase in $\kappa$ should result in increased growth including maximum length and decreased reproduction.

default(leg=false)

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:e_50_adpt] = 1e-3
flea_params[:beta_e_adpt] = 10.
flea_params[:TKTD][CU].TK[KP].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
    LL3GM,
    (100., 2., 1.)
)

outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450), bottommargin=5mm)

Decrease in $\kappa$ should result in decreased growth and (slightly) increased reproduction.
The stimulating effect on reproduction should be relatively subtle compared to the stimulating effect on growth resulting from increased $\kappa$, because it implies a gross decrease in growth, therefore a gross decrease in reproduction, as well as a gross increase in reproduction. The parameter combination below should lead to a net decrease in reproduction.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
    LL3GM,
    (100., 2., -.5)
)

outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

In the code below, the TK parameter is lowered which should produce a hormetic effect on reproduction.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
    LL3GM,
    (100., 2., -.5)
)

outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))

outs = LifeTableDataset[]
C_sweep = geomrange(0.1, 1000.)
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

LifeTableDataset(11000×60 DataFrame
   Row │ t_day  unique_id  species_idx  kappa     kappa_EX  q_p_rel  J_X       ⋯
       │ Any    Any        Any          Any       Any       Any      Any       ⋯
───────┼────────────────────────────────────────────────────────────────────────
     1 │ 0.0    1          1            0.66472   1.0       1.0      0.222054  ⋯
     2 │ 1.0    1          1            0.672075  1.0       1.0      0.459321
     3 │ 2.0    1          1            0.668482  1.0       1.0      0.758668
     4 │ 3.0    1          1            0.667277  1.0       1.0      1.10867
     5 │ 4.0    1          1            0.666679  1.0       1.0      1.49666   ⋯
     6 │ 5.0    1          1            0.666323  1.0       1.0      1.91205
     7 │ 6.0    1          1            0.666089  1.0       1.0      2.34598
     8 │ 7.0    1          1            0.665924  1.0       1.0      2.79109
   ⋮   │   ⋮        ⋮           ⋮          ⋮         ⋮         ⋮        ⋮      ⋱
 10994 │ 15.0   1          1            0.33804   1.0       1.0      2.16702   ⋯
 10995 │ 16.0   1          1            0.337747  1.0       1.0      2.24583
 10996 │ 17.0   1          1            0.337501  1.0       1.0      2.31961
 10997 │ 18.0   1          1            0.337291  1.0       1.0      2.3886
 10998 │ 19.0   1          1            0.337111  1.0       1.0      2.45305   ⋯
 10999 │ 20.0   1          1            0.336957  1.0       1.0      2.51317
 11000 │ 21.0   1          1            0.336823  1.0       1.0      2.56923
                                               53 columns and 10985 rows omitted, 11000×5 DataFrame
   Row │ t_day  rep    L          carapace_length  C       
       │ Any    Int64  Any        Float64          Float64 
───────┼───────────────────────────────────────────────────
     1 │ 0.0        1  0.0264072         0.11949       0.1
     2 │ 1.0        1  0.037256          0.168579      0.1
     3 │ 2.0        1  0.047536          0.215095      0.1
     4 │ 3.0        1  0.0572316         0.258967      0.1
     5 │ 4.0        1  0.0663241         0.300109      0.1
     6 │ 5.0        1  0.0748312         0.338603      0.1
     7 │ 6.0        1  0.0827806         0.374573      0.1
     8 │ 7.0        1  0.0902034         0.40816       0.1
   ⋮   │   ⋮      ⋮        ⋮             ⋮            ⋮
 10994 │ 15.0      10  0.0791303         0.358056   1000.0
 10995 │ 16.0      10  0.0805421         0.364444   1000.0
 10996 │ 17.0      10  0.0818418         0.370325   1000.0
 10997 │ 18.0      10  0.0830385         0.37574    1000.0
 10998 │ 19.0      10  0.0841408         0.380727   1000.0
 10999 │ 20.0      10  0.0851562         0.385322   1000.0
 11000 │ 21.0      10  0.0860918         0.389556   1000.0
                                         10985 rows omitted, 1100×4 DataFrame
  Row │ t_day    num_surviving  survival  C       
      │ Float64  Int64          Float64   Float64 
──────┼───────────────────────────────────────────
    1 │     0.0             10       1.0      0.1
    2 │     1.0             10       1.0      0.1
    3 │     2.0             10       1.0      0.1
    4 │     3.0             10       1.0      0.1
    5 │     4.0             10       1.0      0.1
    6 │     5.0             10       1.0      0.1
    7 │     6.0             10       1.0      0.1
    8 │     7.0             10       1.0      0.1
  ⋮   │    ⋮           ⋮           ⋮         ⋮
 1094 │    15.0             10       1.0   1000.0
 1095 │    16.0             10       1.0   1000.0
 1096 │    17.0             10       1.0   1000.0
 1097 │    18.0             10       1.0   1000.0
 1098 │    19.0             10       1.0   1000.0
 1099 │    20.0             10       1.0   1000.0
 1100 │    21.0             10       1.0   1000.0
                                 1085 rows omitted)

@df @subset(out.repro, :t_day.==maximum(:t_day)) lineplot(
    :C, :cum_repro, xscale=:log10, marker=true, xlabel="External concentration", ylabel="21-day reproduction",
    title=L"PMoA $KP$: Effect on 21-day reproduction"
    )

Embryonal investment

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
    LL3AR,
    (100.0, 2.0, 0.5)
)

outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,200))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450), bottommargin=5mm)

# the same for more concentrations
outs = LifeTableDataset[]
C_sweep = geomrange(0.1, 1000.)
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

LifeTableDataset(11000×60 DataFrame
   Row │ t_day  unique_id  species_idx  kappa     kappa_EX  q_p_rel  J_X       ⋯
       │ Any    Any        Any          Any       Any       Any      Any       ⋯
───────┼────────────────────────────────────────────────────────────────────────
     1 │ 0.0    1          1            0.66472   1.0       1.0      0.222054  ⋯
     2 │ 1.0    1          1            0.672075  1.0       1.0      0.459321
     3 │ 2.0    1          1            0.668482  1.0       1.0      0.758668
     4 │ 3.0    1          1            0.667277  1.0       1.0      1.10867
     5 │ 4.0    1          1            0.666679  1.0       1.0      1.49666   ⋯
     6 │ 5.0    1          1            0.666323  1.0       1.0      1.91205
     7 │ 6.0    1          1            0.666089  1.0       1.0      2.34598
     8 │ 7.0    1          1            0.665924  1.0       1.0      2.79109
   ⋮   │   ⋮        ⋮           ⋮          ⋮         ⋮         ⋮        ⋮      ⋱
 10994 │ 15.0   1          1            0.665441  1.0       1.0      6.21733   ⋯
 10995 │ 16.0   1          1            0.665421  1.0       1.0      6.59176
 10996 │ 17.0   1          1            0.665405  1.0       1.0      6.95057
 10997 │ 18.0   1          1            0.665391  1.0       1.0      7.29334
 10998 │ 19.0   1          1            0.665382  1.0       1.0      7.62018   ⋯
 10999 │ 20.0   1          1            0.665374  1.0       1.0      7.9311
 11000 │ 21.0   1          1            0.665369  1.0       1.0      8.22618
                                               53 columns and 10985 rows omitted, 11000×5 DataFrame
   Row │ t_day  rep    L          carapace_length  C       
       │ Any    Int64  Any        Float64          Float64 
───────┼───────────────────────────────────────────────────
     1 │ 0.0        1  0.0264072         0.11949       0.1
     2 │ 1.0        1  0.037256          0.168579      0.1
     3 │ 2.0        1  0.047536          0.215095      0.1
     4 │ 3.0        1  0.0572316         0.258967      0.1
     5 │ 4.0        1  0.0663241         0.300109      0.1
     6 │ 5.0        1  0.0748312         0.338603      0.1
     7 │ 6.0        1  0.0827806         0.374573      0.1
     8 │ 7.0        1  0.0902034         0.40816       0.1
   ⋮   │   ⋮      ⋮        ⋮             ⋮            ⋮
 10994 │ 15.0      10  0.134172          0.607113   1000.0
 10995 │ 16.0      10  0.138129          0.625018   1000.0
 10996 │ 17.0      10  0.141816          0.641704   1000.0
 10997 │ 18.0      10  0.145253          0.657253   1000.0
 10998 │ 19.0      10  0.148455          0.671742   1000.0
 10999 │ 20.0      10  0.151439          0.685242   1000.0
 11000 │ 21.0      10  0.154218          0.697821   1000.0
                                         10985 rows omitted, 1100×4 DataFrame
  Row │ t_day    num_surviving  survival  C       
      │ Float64  Int64          Float64   Float64 
──────┼───────────────────────────────────────────
    1 │     0.0             10       1.0      0.1
    2 │     1.0             10       1.0      0.1
    3 │     2.0             10       1.0      0.1
    4 │     3.0             10       1.0      0.1
    5 │     4.0             10       1.0      0.1
    6 │     5.0             10       1.0      0.1
    7 │     6.0             10       1.0      0.1
    8 │     7.0             10       1.0      0.1
  ⋮   │    ⋮           ⋮           ⋮         ⋮
 1094 │    15.0             10       1.0   1000.0
 1095 │    16.0             10       1.0   1000.0
 1096 │    17.0             10       1.0   1000.0
 1097 │    18.0             10       1.0   1000.0
 1098 │    19.0             10       1.0   1000.0
 1099 │    20.0             10       1.0   1000.0
 1100 │    21.0             10       1.0   1000.0
                                 1085 rows omitted)

PMoA Embryonal investement leads to an increase in cumulative reproduction.

p = @df @subset(out.repro, :t_day.==maximum(:t_day)) plot(
    lineplot(
        :C, :E_o, xscale=:log10, xlabel="External concentration", ylabel="Energy investment per egg (J)", lw=2
    ),
    lineplot(
        :C, :cum_repro ./ 100., xscale=:log10,  xlabel="External concentration", ylabel="Control-normalized \n 21-day reproduction",
        lw=2,
        )
)
savefig(plot(p, dpi=300), "figures/pmoa_E_repro.png")
p

Each embryo starts with less energy, so there should be implications for neonate fitness.

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [100.]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
    LL3AR,
    (100.0, 2.0, 0.5)
)

DRCModel(LL3AR, (100.0, 2.0, 0.5))

global_params[:remove_juveniles] = []

Any[]

neonate_fitness = DataFrame()

C_sweep = geomrange(0.1, 1000., length=10)
for C in C_sweep
    global_params[:C][1] = C
    m = initialize_model(global_params, phyto_params, [flea_params])
    run!(m)
    all_daphnids = daphnid_output(m, Inf)
    all_daphnids[:,:C] .= C
    append!(neonate_fitness, all_daphnids)
end

Different state variables at birth are plotted below.

PMoA embryonic investment decreases energy at birth and length at birth.
Scaled reserves are barely affected, and age at birth is slightly decreased ()

df = combine(groupby(@subset(neonate_fitness, :unique_id.>1), [:unique_id, :C])) do df
    birth = @subset(df, :t_day.==minimum(:t_day))
    DataFrame(
        e_b = birth.e,
        E_b = birth.E,
        L_b = birth.L,
        a_b = birth.age
    )
end 

 @df df plot(
    lineplot(:C, :E_b, ylim=(0, maximum(:E_b)), xscale=:log10, ylabel="Reserves at birth"),
    lineplot(:C, :e_b, ylim=(0, maximum(:e_b)), xscale=:log10, ylabel="Scaled reserves at birth"),
    lineplot(:C, :L_b, ylim=(0, maximum(:L_b)), xscale=:log10, ylabel="Length at birth"),
    lineplot(:C, :a_b, ylim=(0, maximum(:a_b)), xscale=:log10, ylabel="Age at birth"), 
    xlabel="Exposure concentration", size=(600,600)
)

Survival

include("src/DEB.jl")

apply_effects!

flea_params[:TKTD][CU].TD

7-element Vector{TDModel}:
 TDModel(DRCModel(LL3GM, (0.0, 0.0, 0.0)), :s_G_z, prod)
 TDModel(DRCModel(LL3GM, (0.0, 0.0, 0.0)), :s_M_z, prod)
 TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_A_z, prod)
 TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_R_z, prod)
 TDModel(DRCModel(LL3GUTS, (0.0, 0.0, 0.0)), :h_z, sum)
 TDModel(DRCModel(LL3AR, (100.0, 2.0, 0.5)), :s_E_o_z, prod)
 TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_kappa_z, prod)

global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]

flea_params = deepcopy(default_fleas.DM_AMP) 
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1

flea_params[:TKTD][CU].TK[H].params = [1.]
flea_params[:TKTD][CU].link_TK = false
flea_params[:TKTD][CU].TD[H].DRC = DRCModel(
    LL3GUTS,
    (100.0, 2.0, 1.0)
)

outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100., 1000.]
for C in C_sweep
    global_params[:C][1] = C
    out = life_table(global_params, phyto_params, flea_params; n_reps=10)
    out.repro.C .= C
    out.growth.C .= C
    out.survival.C .= C
    push!(outs, out)
end
out = concat(outs)

sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
    :t_day, :carapace_length, group=:C, 
    xlim=(-1,22), ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...), 
    size=(800,300), lw=1.5
)
pr = @df out.repro plot(
    :t_day, :cum_repro, group=:C, 
    layout=(1,4), xlim=(-1,22), ylim=(0,150)
    )
ps = @df out.survival plot(
    :t_day, :survival, group=:C, 
    layout=(1,4), ylim=(0,1.01)
    )
plot(
    pg, pr, ps, 
    layout=(3,1), xticks=0:7:21, 
    xlabel="Time (d)", lw=2, xlim=(-1,22),
    size=(1000, 450), bottommargin=5mm
    )